#include <cstdlib>

#include <vector>
#include <iostream>
#include <fstream>
#include <string>
#include <algorithm>
#include <boost/timer/timer.hpp>
#include <boost/program_options.hpp>

namespace po = boost::program_options;
using std::vector;
using std::string;

po::options_description desc("Options");

void die(const char *msg) {
	std::cout << msg << std::endl;
	desc.print(std::cout);	
	exit(EXIT_FAILURE);
}

enum class BasePair {
	A, T, C, G
};

std::ostream& operator<<(std::ostream& os, const BasePair& pair) {
	switch(pair) {
	case BasePair::A:
		os << "A"; break;
	case BasePair::T:
		os << "T"; break;
	case BasePair::C:
		os << "C"; break;
	case BasePair::G:
		os << "G"; break;
	}
	return os;
}

int hamming_distance(vector<BasePair>::const_iterator dna_it, 
					 vector<BasePair>::const_iterator lmer_start, 
					 vector<BasePair>::const_iterator lmer_end) {
	int distance = 0;
	std::for_each(lmer_start, lmer_end, [&distance, &dna_it](const BasePair val) {
		distance += *dna_it++ != val;
	});
	return distance;
}

int total_distance(const vector<vector<BasePair> >& dna, 
					 vector<BasePair>::const_iterator lmer_start, 
					 vector<BasePair>::const_iterator lmer_end, 
					 int lmer_size) {
	int total_distance = 0;
	std::for_each(dna.begin(), dna.end(), 
			[&lmer_start, &lmer_end, lmer_size, &total_distance]
			(const vector<BasePair>& sequence) {
		int min_distance = INT_MAX;
		auto it = sequence.begin();
		// position after last possible starting position
		auto end = sequence.end() - lmer_size + 1;
		while (it != end) {
		      min_distance = std::min(min_distance, hamming_distance(it++, lmer_start, lmer_end));
		}
		total_distance += min_distance;
	});
	return total_distance;
}

bool bypass_vertex(vector<BasePair>& sequence, int& pos) {
	while (pos >= 0 && sequence[pos] == BasePair::G) 
		pos--;
	if (pos == -1) return false;
	sequence[pos] = (BasePair) ((int)sequence[pos] + 1);
	return true;
}

bool next_sequence(vector<BasePair>& sequence, int& pos) {
	if ((unsigned)pos == sequence.size() - 1) {
		// leaf node - "bypass" to next node
	        return bypass_vertex(sequence, pos);
	} else {
		// go down further tree
		pos++;
		sequence[pos] = BasePair::A;
	}
	return true;
}

const vector<std::pair<vector<BasePair>,int> > find_motifs(const vector<vector<BasePair> >& dna, int lmer_length) {
	// position of last set value
	int pos = 0;
	vector<BasePair> sequence(lmer_length);
	int best_distance = INT_MAX;
	vector<std::pair<vector<BasePair>,int> > best_words;
	while(true) {
		int distance = total_distance(dna, sequence.begin(), sequence.begin() + pos + 1, lmer_length);		
		if (pos == lmer_length - 1) {
			// leaf node
			if (distance <= best_distance) {
				if (distance < best_distance) {
					best_distance = distance;
					best_words.clear();
				}
				std::pair<vector<BasePair>,int> AnnotatedSequence(sequence,distance);
				best_words.push_back(AnnotatedSequence);
			}
			if (!next_sequence(sequence, pos)) break;
		} else {
			if (distance > best_distance) {
				if (!bypass_vertex(sequence, pos)) break;
			} else {
				if (!next_sequence(sequence, pos)) break;
			}
		}
	}
	return best_words;
}

const vector<vector<BasePair> > parse_file(const string& filename) {
	std::ifstream file(filename, std::ios::in);
	if (file.fail()) die("Could not open input file.");
	vector<vector<BasePair> > data;
	std::string line;
	bool first = true;
	unsigned sequence_size;
	while (std::getline(file, line)) {
		vector<BasePair> sequence(line.size());
		auto it = sequence.begin();
		std::for_each(line.begin(), line.end(), [&it](const char val) {
			switch(tolower(val)) {
			case 'a':
				*it++ = BasePair::A; break;
			case 't':
				*it++ = BasePair::T; break;
			case 'c':
				*it++ = BasePair::C; break;
			case 'g':
				*it++ = BasePair::G; break;
			default:
				die("Invalid alignment character"); break;
			}
		});
		if (first) {
			if (sequence.size() == 0) die("DNA sequence length cannot be zero.");
			sequence_size = sequence.size();
		} else {
			if (sequence.size() != sequence_size)
				die("All DNA sequences have to have same length.");
		}
		data.push_back(sequence);
	}
	return data;
}

void print_to_file(string filename, const vector<std::pair<vector<BasePair>,int> >& sequences) {
	std::ofstream file(filename, std::ios::trunc | std::ios::out);
	if (file.fail()) die("Could not open output file.");
	std::for_each(sequences.begin(), sequences.end(), [&file](const std::pair<vector<BasePair>,int>& sequence) {
	    std::for_each(sequence.first.begin(), sequence.first.end(), [&file](BasePair pair) {
			file << pair;
		});
                file << " - " << sequence.second;
		file << std::endl;
	});
}

int main(int argc, char **argv) {
	string inputfile, outputfile;
	int motif_length;
	bool print_time;
	desc.add_options()
		("inputfile,I", po::value(&inputfile), 
			"input filename. File has to contain sequences of base pairs without spaces, sequences separated by newline.")
		("motif_length,L", po::value<int>(&motif_length), "Length of best motif to find.")
		("timing,T", po::value<bool>(&print_time)->default_value(false), "Print timing (user + system) in ms to stdout.")
		("outputfile,O", po::value(&outputfile),
			"output filename. File contains all optimal motifs.")
	;	
	po::variables_map vm;
	po::store(po::parse_command_line(argc, argv, desc), vm);
	po::notify(vm);
	if (!vm.count("inputfile")) 
		die("No inputfile specified.");
	if (!vm.count("motif_length")) 
		die("No motif length specified.");
	
	const vector<vector<BasePair> > dna_sequences = parse_file(inputfile);
	boost::timer::cpu_timer timer;
	const vector<std::pair<vector<BasePair>,int> > motifs = find_motifs(dna_sequences, motif_length);
	if (print_time) {
		auto time = timer.elapsed();
		// we want time in ms
		std::cout << (time.system + time.user) / 1000000 << std::endl;
	}
	if (vm.count("outputfile") != 0) {
		print_to_file(outputfile, motifs);
	}
	return EXIT_SUCCESS;
}


